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Models that invoke nonlinear wavefront propagation in a chemically excitable medium are rife 

^4— I ■ in the biological literature. Indeed, the idea that wavefront propagation can serve as a signaling 

O , 

. , mechanism has often been invoked to explain synchronization of developmental processes. In this 



paper we suggest a new kind of signaling based not on diffusion of a chemical species but on the 
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T^ • propagation of mechanical stress. We construct a theoretical approach to describe mechanical 
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Q . signaling as a nonlinear wavefront propagation problem and study its dependence on key variables 



such as the effective elasticity and damping of the medium. 



[C . PACS numbers: 87.17.Aa, 87.53.Ay 

in 

en 

^' 

o 



X 



The physical phenomenon of nonhnear wavefront propagation in an excitable medium is 
widely exploited by biological systems to transmit signals across many cells. For example, 
when the slime mold Dictyostelium begins to aggregate to form a fruiting body, wavefronts 
of the molecule cAMP propagate across the amoeba colony [1-3]. Although cAMP itself 
spreads diffusively, wavefronts of cAMP propagate ballistically across the amoeba colony 
because the colony is chemically excitable: when the local concentration of cAMP exceeds a 
threshold, further local release of the species is triggered [2, 3]. Similar wavefronts of calcium 
and potassium, respectively, signal fertilization in eggs [4] and the onset of spreading cortical 
depression [5], associated with migraine auras. 

In recent years, however, there has been a growing recognition that mechanics plays an 
important role in biology, and that many cells sense and respond not only to chemical stimuli 
but also to mechanical stimuli [6-10]. This raises the possibility that mechanosensing at the 
cellular level could give rise to collective phenomena at larger length scales such as collective 
cell migration [11, 12]. 

In this paper, we consider the possibility of mechanically-induced waves at the tissue level. 
Biological systems are typically overdamped, so they do not support sound waves and stress 
cannot propagate ballistically. However, mechanical signaling via ballistic propagation of a 
nonlinear wavefront can occur in a mechanically-excitable medium, much as chemical sig- 
naling via ballistic propagation of a nonlinear wavefront can occur in a chemically-excitable 
medium. It has recently been suggested that two different biological systems might be 
mechanically-excitable: the early Drosophila erahrjo [13] and the developing heart [14]. The 
early Drosophila embryo supports mitotic wavefronts: nuclei at the poles of the embryo tend 
to divide first, giving rise to a mitotic wavefront separating dividing nuclei from those that 
have not yet divided. This wavefront propagates across the entire embryo [15]. Likewise, 
the heart tube of the avian embryo beats via contractile wavefronts that are initiated at 
one end of the tissue and propagate across the heart tube with each beat [16]. These two 
examples are very different in biological details but share a key common feature-that the 
nuclei in the case of the Drosophila embryo and cardiomyocytes in the case of the heart 
tube both generate stresses as they proceed through mitosis and contraction, respectively. 
It is therefore important to develop a general physical understanding of the basic properties 
needed to yield mechanically-induced nonlinear wavefronts. 

Here we introduce two minimal models of mechanically-excitable media. In each case. 



stress can be released at certain sites, or nodes, if the local stress exceeds some threshold 
value. In the case of the early Drosophila embryo, these nodes would represent cell nuclei, 
while in the case of the developing heart tube, the nodes would represent cardiomyocytes 
(heart cells that can contract). In our models, the released stress is transmitted through 
the damped elastic medium, potentially causing further release of stress at other nodes. We 
solve these models and identify characteristic features exhibited by nonlinear wavefronts in 
such systems. 

We start by considering two simple examples of media that can support mechanical stress 
and have overdamped mechanics. The elasticity is characterized by Lame coefficients A 
and fi, or equivalently by the Young's modulus E and dimensionless Poisson ratio u within 
linear elasticity theory [17]. These parameters relate the stress aij inside the elastic material 
to its strain (deformation) eij: 

E r V 1 

where d is the number of dimensions (2 or 3), and summation over repeated indices is 
implied. The strain Sij is defined in therms of the displacement vector Ui of the elastic 
material: Eij = ^{djUi + diUj). To avoid confusion we will follow the usual convention and 
label the two-dimensional versions of the parameters with a subscript 2, they are related to 
their three-dimensional counterparts by E2 = E/{1 — z/^) and z/2 = z//(l — z/) [17]. The force 
per unit area is given by the divergence of the stress: Pi = djCJij. 

In the simplest model that we consider, corresponding to a thin elastic film that slides 
frictionally over a surface, we balance this elastic force with a friction term, TdtUi, where F 
is the friction coefficient. Such a system can be described by a two-dimensional model with 
the equation of motion: 



TdtUi 



djdjUi + didjUj 
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Note that Eq. 2, which describes the response of elastic medium to a displacement Ui, is 
similar to the diffusion equation, but is a tensor equation instead of a scalar one. 

The second model that we consider is a three-dimensional realization of an overdamped 
elastic medium, such as a polymer network immersed in a fluid. The elasticity of the system 
is also described by Eq. 1, but the friction force is proportional to the relative motion of the 
fluid and the elastic network: T{dtUi — Vi). The stress in an incompressible fluid depends 



linearly on the pressure p and the shear rate, 7™*^ = \{diVj + djVi): 

^r = -P^^3 + 2r/7r^ (3) 

in both two and three dimensions [18]. In the over damped limit (zero Reynolds number), 
taking the divergence of (3) gives the Stokes equation. Combining the elastic and fluid 
equations gives a closed system for -Uj, Vi and p: 

E 

(rf-l)z/ 



T{dtUi - Vi 



djdjUi + — — T^^didkUk 



(4) 



2(1 + z/ 
T{dtUi - Vi) = diP - r]djdjVi, (5) 
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Eqs. 4-6 are identical to the two-fluid model studied by Levine and Lubensky [19], but 
without the inertial terms. 

We now add mechanical excitability as follows. We consider a collection of nodes at 
positions {Rn}, where n indexes the nodes. A node can be activated if some measure of the 
stress (for example, the absolute value of its largest eigenvalue) exceeds a threshold value a. 
If this occurs at time t, the node releases additional stress over a time interval At. For a 
node at Rn activated at time t = tn, we therefore introduce an extra force into equation (2), 
of the form 

^.active _ Q.Q^.S^^ _ ^^)0(t _ t„)0(t„ + At - t), (7) 

where Qij is a tensor of rank 2, which in general can have three contributions: a hydrostatic 
expansion/contraction, a rotation, and a traceless dipole. 

We solve for the response of the two-dimensional overdamped elastic medium of Eq. 2 
to the active force in Eq. 7 by deriving the Green's tensor Gijk{x,t), which relates the 
displacement Uk{x,t) to a source term Qij6{x)Q(t) at the origin at time t = 0. We find 
that the material parameters E2, V2 and F combine in two quantities with the dimensions 
of diffusion constants, 
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which correspond to motion in the longitudinal and transverse directions respectively, and 
together completely determine the solution. Here /x = E2/2{1 + V2) = E/2{1 + u) is the 
material's shear modulus, which is the same in two and three dimensions. The resulting 



Green's tensor is given by: 



Gijk{x,t) = - 
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where x = \/x ■ x and 
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We can derive a similar solution for the response of the two-fluid model of Eqs. 4-6 to 
the active force in Eq. 7. In this case there is an extra parameter, the viscosity rj of the 
fluid, which gives rise to a natural relaxation timescale r = n/rj of the system. In three 
dimensions, the three quantities governing the solution of equations (4-6) are given by: 
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For the associated Green's tensor we find: 

(j'ijkyx, t) = Gjjfc [x, t) + (jijk \x, t), 
1 
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For an input term that runs only over a time interval At as in Eq. 7, continuity demands 
that for t < At we have Gijk{x,t) = G^°^{x,t) + Gl^f{x,t), and for t > At this changes to 
G,,fc(f , t) = G'l^^ix, t) - G^^^ix, t - At). 



3y I Trerf ( 

2 V V2v^ 



V^/^^- Trerf ^ ^ 



Because our model equations are linear, we can now use the principle of superposition to 
study the effect of many source terms. We initialize the system by activating a single node at 
the origin at t = 0. We then measure the stress at the other nodes as a function of time, and 
activate them if they are above threshold by releasing more stress, according to equation (7). 
We consider various cases for the arrangement of the nodes: a regular triangular lattice, a 
random configuration with short-range correlations (as in a random packing of disks) and an 
uncorrelated random configuration. In addition, we look at variants in which the force term is 
purely isotropic (hydrostatic expansion/contraction) or is in the form of a volume-conserving 
force dipole, with either random orientation or orientations correlated to the direction of 
the traveling wavefront. In all cases, the model produces an activation wavefront with a 
well-defined speed, as shown in Fig. 1. We find that the speed of the wavefront depends 
on the density of nodes but is insensitive to their arrangement. However, the spread of 
the wavefront around its mean increases with the amount of randomness (Fig. 1(b)). Not 
surprisingly, if the orientations of the force dipoles are chosen at random, the speed of the 
wavefront is the same in all directions so that its shape is circular, as in Fig. l(b, inset). 
In contrast, if the dipoles are all oriented in the same direction, the wavefront is no longer 
uniform, but is faster in the direction of orientation. Also, for the same magnitude of 
the active force, the speed of the wavefront is somewhat higher if the force dipole Qij is 
hydrostatic than if it is a traceless dipole (Fig. 1(b)). All these observations indicate that 
the wavefront speed is dictated primarily by the properties of the medium and the average 
distance between nodes, and is insensitive to both the form of the active force and spatial 
arrangement of nodes. 

Dimensional analysis of the material parameters of our two-dimensional model (2) shows 
that there is only one possible scaling for the wavefront speed with the material parameters: 
V ~ E2/{aV), where a is the grid spacing. The dimensionless speed v = {aV / E2)v depends 
on the material's Poisson ratio V2 and the dimensionless threshold a = a* a? /Q^ where Q 
is the strength of the force term. We have determined the function ■Z;(z/2,q; numerically for 
both hydrostatic and dipole force terms. We find that it obeys a fairly simple functional 
form, which can be motivated by an analytical argument based on the case of the simplest 
force term, the purely hydrostatic one Qij = Q6ij. In this case, the stress is given by 
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FIG. 1. Calculated wavefronts. (a) Plots showing for each node (dots) the distance from the first 
activated node vs. time of activation, with linear fits. This example has a dipole source term 
and shows results for three different types of grids, all with the same density: regular triangular 
(gold), correlated random (blue) and uncorrelated random (red), (b) Mean wavefront speed for 
six realizations, with hydrostatic force term (1-3, blue) and dipole force terms (4-6, red), and on 
three different grids: regular triangular (circle, 1 & 4), correlated random (diamond, 2 &; 5), and 
uncorrelated random (square, 3 &: 6). Error bars indicate standard deviations. Inset shows the 2- 
dimensional field of nodes, indicating the orientation in which nodes are activated, and color-coded 
according to the time at which they are activated, on a hue scale (red-yellow-green-blue-violet). 

Because the stress drops off quadratically with distance, the major contribution to the stress 
at any node is due to forces exerted by neighboring nodes. Moreover, since the front expands 
radially, typically only a single nearest neighbor of any node will have been activated recently. 
We can therefore get a reasonable estimate for the local stress at a node by considering that 
nearest neighbor to be the only source. We introduce the dimensionless time i = E2t/{a'^r)] 
then for a single source a distance a away, the time at which the largest eigenvalue of the 
stress (19) reaches the dimensionless threshold a is given by: 

1 + 1^2 



a 



1 - ^2 



2t 



,-(l-z.|)/4t 



(20) 



Unfortunately, Eq. (20) cannot be inverted analytically. However, the two factors containing 
i are easily inverted, allowing us to make an educated guess for the functional form of the 
resulting dimensionless speed: 

4(ciQ; + C2) log(Q;) 
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where c\ and C2 need to be determined numerically; we find c\ = 4.0 and C2 = 1.5. As 
shown in Fig. 2(a), the form given by equation (21) works remarkably well. Moreover, the 
same functional form also describes the results for a dipole force term wavefront, as shown 
in Fig. 2(b), the only difference being the values of the two fit parameters - here we find 
Ci = —1.0 and C2 = 1.0. 

In line with intuition, our model predicts that there is a maximum threshold value cimax 
above which a wavefront will not propagate. This can happen for one of two reasons: either 
the force is not large enough to create a stress at the next node that exceeds the threshold 
value, or the nodes are so far apart that, due to the diffusive nature of the stress spreading, 
the threshold value is not reached. Both possibilities are contained in the form of the 
dimensionless version of cimax, given by the maximum of the right hand side of equation (20), 
which gives cimax = 2e^^-^^'^^^^'^ at tmax = ^"Y^ and corresponding to a minimum speed 



"^min = i~[r- ^^ note that as z/2 approaches its maximum value of 1, the minimum speed 
diverges, as can be seen in Fig. 2(a). 

For the three-dimensional two-fiuid model of Eqs. (4-6), there are two independent quanti- 
ties with the dimensions of speed, E/aT and a/r, where r = rj/fi is the material's relaxation 
time (Eq. 12). We note that both of these scale linearly with the material's Young's modu- 
lus E (or equivalently, with the material's shear modulus ^), which implies that also in this 
case the resulting wavefront velocity in a similar setup with excitable nodes will scale linearly 
with that modulus. It will also scale with F^^r^^^", where n is some number between and 
1, indicating that both the internal viscosity of the moving fiuid and the friction between 
the elastic and viscous material contribute to the damping of the ballistic motion. 

In this paper, we have introduced two theoretical realizations of mechanical signaling 
in an overdamped elastic medium. We have shown that nonlinear wavefront propagation 
in the models is a robust feature of both models. In both cases, the wavefront velocity is 
insensitive to the spatial distribution of excitable nodes. It is also insensitive to whether 
the stress is released in an isotropic or traceless anisotropic fashion. Furthermore, a funda- 
mental feature of both models is that the wavefront velocity is proportional to the Young's 
modulus of the medium, and the magnitude of the velocity can be understood simply and 
quantitatively in terms of characteristic dimensionless variables, such as the stress threshold 
made dimensionless with the magnitude of the force dipole released when a node is excited. 

The overdamped elastic models considered here are the simplest models that could be used 
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FIG. 2. Dimensionless wavefront speed as a function of the Poisson ratio V2 and dimensionless 
threshold a. Symbols indicate numerical solutions of the full system, lines the functional form 
of equation (21). (a) Hydrostatic force term Qij = Q5ij. Fit parameters ci = 4.0, C2 = 1.5. 
Values of a: 0.1 (blue/dots), 0.2 (red/ diamonds), 0.3 (gold/squares) and 0.4 (green/triangles), 
(b) Dipole force term with random orientation angle 9: Qij = —Q cos{29){6ii6ji — 8128 j2) — 
Qsm.{26)[5ii5j2 + 5i25ji). Fit parameters ci = —1.0, C2 = 1.0. Values of a: 0.050 (blue/dots), 
0.075 (red/diamonds), 0.100 (gold/squares), 0.125 (green/triangles up), 0.150 (purple/triangles 
down), and 0.200 (pink/hexagons). 



to describe a tissue. It would be worthwhile to explore mechanical signaling in other models 
that have been proposed for tissues, including the active gel model [11, 20, 21] and cellular 
models [22, 23]. The active gel model of Kruse et al. [20, 21] is an extension of the two-fluid 
model of Levine and Lubensky [19] with a continually active (energy-consuming) term to 
model the dynamics of the cytoskeleton due to motor activity. As we have shown here, such a 
continuous activity is not necessary to describe wavefront propagation, as local and discrete 
activity is sufficient. However, given the presence of active motors in the cytoskeleton, it 
would be interesting to see how the wavefront is affected by an active term in the model. It 
would also be interesting to compare the results of such an interaction with those of Bois 
et al. who study pattern formation in active fluids due to chemical signaling [11]. These 
active models, and the one we used here, are continuous models. However, tissues are of 
course composed of discrete units, the cells. As shown by Manning et al. several mechanical 
properties of the tissue, like its surface tension, are determined by cell-cell adhesion and 



cortical tension [22] . Recent work by Chiou et al. provides a method to measure the relative 
magnitude of forces acting within and between cells [23] . These results now make it possible 
to construct a quantitative cell-based tissue model, in which wavefront propagation due to 
mechanical signaling can be studied as well. 

Now that we have introduced a minimal model for mechanical signaling via nonlinear 
wavefront propagation, we can ask how one might identify biological contexts in which me- 
chanical signaling is likely to occur. We can also ask how to determine whether a given 
wavefront is an example of mechanical signaling. Wavefronts of processes that generate 
stresses are obvious likely candidates. In order for a medium to be mechanically-excitable, 
however, it is not enough to have a collection of nodes capable of generating stress. There 
must also be a mechanism to trigger the nodes to release stress once a stress threshold is 
reached. A likely mechanism would involve stress-dependent ion channels that can release 
ions above a threshold stress [24-27]. Thus, a system with known wavefronts of ion con- 
centration would be a possible candidate, especially if the wavefronts can also be triggered 
mechanically. Mitotic wavefronts in the early Drosophila embryo fit many elements of this 
profile [13]. The process of mitosis generates stresses as chromosomes condense and as they 
segregate [28]. There is a known calcium wavefront that propagates across the embryo in 
tandem with the mitotic wavefront [29] . Similarly, contractile wavefronts in heart tissue may 
be a form of mechanical signaling [14]; they generate stresses as cardiomyocytes contract 
in a process that is well known to involve calcium via the excitation-contraction mecha- 
nism [30, 31]. Another possible realization is spreading cortical depression, which involves 
a potassium wavefront [5] and that can be triggered mechanically [32]. These examples 
suggest that it may be worthwhile to re-examine other known examples of ion signaling in 
biological contexts to see if they are more properly interpreted as mechanical signaling. 
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